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Abstract 

The GOY model is a model for turbulence in which two conserved quan- 
tities cascade up and down a linear array of shells. When the viscosity pa- 
rameter, z/, is small the model has a qualitative behavior which is similar to 
the Kolmogorov theories of turbulence. Here a static solution to the model 
is examined, and a linear stability analysis is performed to obtain response 
eigenvalues and eigenfunctions. Both the static behavior and the linear re- 
sponse show an inertial range with a relatively simple scaling structure. Our 
main results are: (i) The response frequencies cover a wide range of scales, 
with ratios which can be understood in terms of the frequency scaling proper- 
ties of the model, (ii) Even small viscosities play a crucial role in determining 
the model's eigenvalue spectrum, (iii) As a parameter within the model is 
varied, it shows a "phase transition" in which there is an abrupt change in 
many eigenvalues from stable to unstable values, (iv) The abrupt change is 
determined by the model's conservation laws and symmetries. 

This work is thus intended to add to our knowledge of the linear response 
of a stiff dynamical systems and at the same time to help illuminate scaling 
within a class of turbulence models. 



PACS: 47.27.-i, 47.27.Jv, 24.27.Eq, 05.45. +b, 02.10 



1 



1 Introduction 



In turbulent flow a hydrodynamic system couples together many different length 
scales and thus shows in a single process a huge range of relaxation rates. There 
are a variety of simplified models of turbulence which are also intended to show this 
wide range of frequency and wave-number scales. One such model goes under the 
inelegant acronym of "GOY". The model couples together a large number of shells, 
each with its characteristic scale of wave-vectors and relaxation times. Shells are 
spaced logarithmically in wave vector. The nth shell is characterized by a single 
complex velocity, U n , which then depends upon time, t. The model is a linked set of 
ordinary differential equations for all these U n with equations picked to mimic those 
of real hydrodynamic flow. 

We do not know whether the model has much to do with turbulence. But it 
certainly illustrates the behavior of a stiff system. (Stiff systems are ones in which 
numerical simulations are made difficult by the effects of a huge range of relaxation 
rates.) In this paper, we describe the time dependence of the model in the simplest 
possible situation, the linear response to disturbances around a static solution. The 
response is described as an eigenvalue problem, with the response matrix being a 
large matrix which inherits the conservation laws, the scaling properties, and the 
symmetry principles of the GOY model. We look for the eigenvalues and eigenstates 
of the matrix. As we find them we see that the linear response, in turn, shows a 
considerable richness including scaling behavior and the analog of a phase transition. 
Much of this behavior can be understood in terms of the several symmetries and 
conservation laws of the original model. 

The richness of behavior in this linear response theory serves to remind us that 
there is one area of scaling or similarity theory which has not been fully explored, 
the determination of eigenf unctions for matrices which have an underlying scaling 
structure. We also do not understand very much about turbulence, or even about 
the flow of information up and down dynamical linear chains. This paper is about 
these three not-fully-understood areas of applied mathematical science. 

To start out, we show the most interesting results of our study. We plot in figure 
|l|, the eigenvalues of the linear stability in a sort of polar diagram in which the polar 
coordinates, 9 and r respectively are the eigenvalues' phase and are proportional to 
the logarithm of the magnitude of the eigenvalue. (See equation fl25|) below for a 
precise definition. )[] The two different parts of the plot show the response to a purely 
real disturbance (in part |l|a) and to a purely imaginary disturbance in part |l]b. This 
distinction is meaningful because both the basic equations and the static solution 
are purely real. For this figure we have picked a particularly small value of the 
viscosity parameter so as to arrive at a simple scaling behavior. Indeed the simple 
spacing of the points in figure la shows that a simple multiplicative law generates 
the higher order eigenvalues from the lower order ones. That is why the points fall 

1 Similar plots are given for different values of e in figures @) and @. Their discussion will be 
postponed to later on in the paper. 
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on a straight line with regular spacings. There is a tremendous range of scaling of 
the eigenvalues, and it all looks provocatively simple. This paper is mostly aimed 
at producing a partial explanation of these figures. 

The paper starts with the introduction of the GOY model, and treats scaling for 
the static solution in section 2. In section 3 we discuss the linear stability analysis in 
general and apply it to the consideration of the scaling of the simpler, real component 
of the response. Also a relation between eigenvalues and conservation laws is derived. 
Section 4 discusses the linear stability in the imaginary response sector and the phase 
transition and scalings seen there. Section 5 summarizes the results. Footnotes and 
appendices discuss findings that are not necessary for the main line of thought. 

1.1 The model 

The basic structure of GOY originated from Gledzer [[]]], and was motivated by the 
cascade structure of turbulent eddies and conservation laws. Later, Ohkitani and 
Yamada 0] generalized the model and carried out numerical studies that revealed 
chaotic behavior and a dynamic scaling of velocity fluctuations. These studies have 
been extended in[^, ||, |5|, |5). A popular introduction can be found in J7J. 

The basic ingredient is the hierarchical structure: The n th shell is characterized 
by a wave vector of length 

k n = k X n , n=l,2,...N, (1) 

with A > 1 and by a complex velocity mode U n . The Navier-Stokes dynamics is 
mimicked by the following set of ODEs: 

^U n = F n -C n (U*)-D n , (2) 
at 

where the three terms represent respectively forcing and cascade processes and dis- 
sipation. (The star * indicates complex conjugation.) 
We pick a forcing on the first shell 

F n = 5 nA f. (3) 

Most previous studies |2|, ||, || El °f t ne GOY model use a forcing on the fourth 
shell. Our choice of the first shell seems to give a simpler structure to the results. 
More details can be found in Appendix A. In our numerical work we shall choose 
A = 2, k = A" 1 , / = 1 unless otherwise stated. 
The dissipation term is 

-D n = -vk 2 n U n (4) 

and is the fc-space representation of the usual viscous dissipation process. 

The cascade term couples the shell n to its nearest and next nearest neighbors, 

C n {U) = k n U n+ iU n+ 2 — ek n _iU n -iU n+ i — (1 — e)k n _2U n -\U n -2- (5) 
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The boundary conditions are simply Z7_i = Uq — and that the velocities go to 
zero as n becomes large. (In our numerics we represent that by cutting off at some 
large shell numbered N and then using the conditions Un+i = Un+2 = 0.) Because 
the static equations couple n with the four immediately neighboring shells, we need 
four boundary conditions to define the problem. Two of the conditions are at the 
low-n end, and two at the high. The nature of the boundary conditions will become 
crucial later on. 

The model parameter e determines the ratio between upscale and downscale 
coupling. It gives us a convenient tool for varying the model and seeing qualitatively 
different ranges of behavior. 

In the inviscid, unforced limit there are two conserved quantities: the energy 

Z n 

and a second conserved quantity || of the form 

\U I 2 

tf = V-^4_ (6b) 

that can be roughly associated with the helicity in fluid motion ||. Even though 
the association is not perfect, we shall call this quantity the helicity in this paper 
I- 

The GOY model shows dynamical as well as static behavior. There is a static 
solution |§ of the GOY model (dll/dt = 0) in which the phases are zero. When we 
wish to point particularly to the static solution, we shall write it as u n , while we 
use U to refer to either the static or the dynamic case. In this paper, we study the 
static solution, u n , and its linear stability properties. 

The phases are not uniquely fixed for u n || which is part of a one parameter 
family of solutions. (We pick the particular solution which has u n real and positive.) 
The invariance in the phase will be important in our linear stability analysis. It gives 
rise to a zero eigenvalue. There are other approximate invariance properties which 
will be discussed when we get to the linear stability analysis. 



2 Scaling (of static solution) 

We are interested in understanding the behavior of the model in the limit as the 
viscosity becomes small. In that limit, there are three regions of fc-space, or of n, 
called the stirring subrange SSR, the inertial subrange ISR, and the viscous subrange 
VSR. The last is dominated by the viscosity term and has a velocity which decays 
very rapidly with k. (In the GOY model, the decay is exponential in a power || of 
k n .) The SSR is naturally enough the range in which stirring is directly important. 
In our case this comprises only n — 1. The inertial subrange is the intermediate 
range of wave vectors between these two. Here, the behavior is dominated by the 
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cascade term. In the center of the ISR, the solution is best described by using the 
product of three velocities: 



A n (u) = k n U n U n+1 U n+2 . (7) 

In the region in which only the cascade term matters, this product is 

A n = A + B(e - l) n (8) 

where A and B are adjustable constants of integration, representing respectively the 
energy and helicity flux though the inertial range. We shall work with e in the unit 
interval so that, for large n, the A term dominates. In this region: 

U n = Wnk- 1 / 3 (9) 

where W n is a periodic function with period three. This behavior continues as n 
increases until one enters the dissipative range, and U n begins to fall off very rapidly. 

Notice that we have, as expected, four undetermined parameters in the solution. 
In the ISR these parameters are A and B and the two parameters which define the 
period-three oscillations. 

The scaling behavior is illustrated in figure [| which plots W n = kj/ 3 U n against n 
for a case in which v has the value 10 _3 16~ 20 . Notice how U n shows a simple scaling 
for large n up to n of about 67, and then it nose-dives. The nose-dive occurs as n 
increases above a dissipative threshold, which is achieved when the dissipation term 
and the cascade terms are roughly of equal size. Usually, W n is of the order unity 
in the inertial range so that this dissipative cutoff can be computed as a function of 
a threshold value of wave vector, called hp, which obeys 

kv i/3 ~ v. (10) 

The condition is thus that the viscous effects should dominate for values of n larger 
than this dissipative threshold, that is, 

3 

n > No ~ --log x u (11) 

There is a simple scaling theory which we can apply to this case. Because the 
static solution has an asymptotic period three ||, the system should be almost 
unchanged by the transformation |: 

v _> v > = z//A 4 (12a) 
N -> N' = N + 3 (12b) 



2 The change in N would be unimportant were we really working with very large values of N. 
For numerical convenience we work with N only seven or eight larger than No ■ This change in iV 



in equation (12a) eliminates effects produced by changing the cutoff. Throughout this paper we 



use a prime to denote quantities changed by the transformation of equations (12a) and (121) 
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In particular, as one goes from the primed to the unprimed situation, the static 
solution should be basically the same at both ends, with the large n-end being only 
modified by having smaller values of u. Let v! be the solution for the velocity in the 
situation changed according to equations (|12a] , |12b|) . The two scaling symmetries 
can be written as: 

(S) u' n = u n {l+0{p-) 2 / 3 ) for N D -n»l, (13a) 

(L) A< +3 = Un(l + 0(e- l) n ) for n » 1 . (13b) 

In these equations "S" stands for small k and "L" refers to large k. Note that the 
first of these equations remains valid in the SSR while the latter remains valid in 
the VSR. 



In equations (|13a|J13b| ) we have included estimates of the error resulting from 
the terms we have not taken into account in constructing the scaling. In the small 
fc-range, we do not include viscous effects, so the error estimate is the relative size of 
the viscous term. For large k, we do not include the helicity flux represented by the 
B in equation (|§D so this effect is put into the error. The global error is set by the 
sizes of each of these errors at the far ends of the viscous subrange. One such error 
is the value of v. The other, and usually larger effect, is the effect of the helicity 
flux term, B, at the large-n end of the ISR. This term has the order of magnitude, 

error = 0(1 - e) N ° . (14) 

This scaling error must be at least as large as the maximum of this error and v. 

In terms of W, defined by equation (||) the two scaling equations imply respec- 
tively that W' n = W n and that W' n = W n+ 3. In the center of the ISR both of these 
scaling symmetries are valid. Then W has a period three symmetry: 

(I) W' n = W n = W n+3 for n » 1 and N D -n»l (15) 



"I" refers to an intermediate range which occurs in the middle of the ISR. 

To check our thinking, we calculate the static solution of the GOY model. De- 
viations from scaling can be seen by looking at the behavior of 

5 s ,n = 1 - (16a) 



W, 



(16b) 

VV n+3 



OL,n 

Both quantities are plotted in figure ||. Also shown in this figure are theoretical 



lines which show the errors defined in equations ( |13a| ),( |I"3b"|) and (14). Theory and 
experiment show excellent agreement. 
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3 Eigenvalue spectrum 



3.1 Linear Response 

The next stage is to do linear stability analysis. We consider small deviations about 
the static solution u n by writing: 

U n {t) = u n {l + 5<b n e at ) (17) 

Since the static solution u is real the eigensolutions split neatly into oscillations 
of the phase and the amplitude, corresponding respectively to 5$ n being real and 
imaginary. Then we have an eigenvalue equation 

a5<5> n = J2 A nm$®m- (18) 

m 

We distinguish the two different cases with subscripts M for magnitude and (f> for 
phase. We use a superscript j to denote which eigenvalue is being considered. In 
what follows we will always order the eigenvalues so that the magnitude of the 
eigenvalue increases with increasing values of the index j. Then the eigenstate will 
be very small for those shells which have n considerably smaller than j. We will 
then argue that the behavior of the response matrix for n and m of order j will play 
a large role in determining the eigenvalue. 

Note that the response matrix A has two parts. The dissipative part is 

D nm = vS nm kl, ( 19a ) 

and the cascade response is: 

d 

C nm = u m - — C n (u). (19b) 
ou m 

Now the structure of A is different for the two kinds of response. For the magnitude 
response, 

A M = -C-D, (20a) 

while for the phase response 

A < f ) = C — D. (20b) 
From (|5p one sees that the matrix — C has rows of the form 

k n (bu n+1 +cu n - 2 )u n -ii 0, k n (u n+2 +bu n -i)u n+1 , k n u n+1 u n+2 , 0, . . . , 

(21) 

where b = — e/A and c = (e — 1)/A 2 . 
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3.2 Eigenvalue spectra 

At the center of the ISR u n shows a simple scaling superposed on top of a period 
three behavior. Thus, the response matrix also has a period three scaling: 

Cn+3,m+3 = A 2 C nim for n, 77i in region I. (22) 

Consequently, if C dominates the behavior of the matrices A in the determination of 
eigenvalues in some range of j, then the eigenvalues would obey the scaling property 

a M 3 = ^ 2<j m f° r J i n region I (23a) 
crj, +3 = A 2 crj, for j in region I (23b) 

In region I, then, the logarithms of the eigenvalues should be evenly spaced along 
straight lines with spacing of log A 2 . Also, we might think that if C dominates the 
behavior of the matrices A, then the phase and magnitude eigenvalues should be 
the same except for a minus sign. (See equations ( |20a| ) and ( |20b| )). The viscosity is 



unimportant in the entire region S. Thus we expect 

a 3 M = — crj, for j in region S. (24) 

Figure [I] shows the eigenvalue spectra for e = 0.3 for both modulus and phase 
stability matrices A M and A^. As the distance between the eigenvalues a- 7 in the 
hierarchical GOY model grows roughly exponentially, it is hard to visualize the 
spectra in the complex plane. To have some kind of visualization we use a kind of 
polar representation in which the phase of the plotted point is exactly the phase of 
cr while the distance from the center of the polar plot is given by: 

r = log A (l + 2 1 V|). (25) 

The factor of 2 10 is put in to enhance the visibility of eigenvalues with small values 
of \a\. The plot has unstable eigenvalues showing up on the right hand side of the 
origin while stable ones show up on the left. For large eigenvalues, even spacings on 
the plot mean that successive eigenvalues have ratios which are a constant. Thus 
even spacings are indicative of some kind of simple scaling. 

In both the magnitude and the phase sector, the eigenvalues are arranged in 
several branches. Within each branch there are regions of even spacing, indicative 
of simple scaling. The magnitude eigenvalues seem particularly simple with three 
well-defined branches: A set of real eigenvalues and a pair of complex conjugate 
branches. All eigenvalues are stable. The phase eigenvalues show a more complex 
structure with what looks like more regions of simple scaling. Nonetheless both 
sets of eigenfunctions show the even spacing demanded by scaling. [] As one might 

3 The existence of three branches is not due to the existence of a period three in the solution. 

— 1/3 

One can assume an approximate solution u n = k n , but the corresponding matrix also consists 
of three radial branches, although we have no period 3 in the "solution" any more. 
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expect, the minimal modulus eigenvalue is about of the size of the entries in the first 
row of the cascade matrix C, the modulus of the second one is of the order of the 
second row, and so on. Generally, the modulus of the nth eigenvalue is about of the 
size of the nth row of the matrix. This rule holds roughly until the eigenvalue with 
j equal to No is reached whereupon the successive eigenvalues are real and have 
values of the order of the diagonal elements in the dissipation matrix, D. In the 
center of the inertial range the eigenvalues change from real to a complex conjugated 
pair and back to real and so on with a period three. 

Thus, much of what we see is what we might expect. But not all. Equation ([23]) 
is completely inconsistent with the pictures we are seeing. This equation implies that 
if we have stable eigenvalues in the magnitude sector, we should have unstable ones 
in the phase sector. But all eigenvalues in figure [I] are stable. This result cannot be 
consistent with the notion that viscosity is unimportant for the eigenvalues in some 
region of the response. We have other difficulties. The pattern of phase response 
eigenvalues looks much more complex than the pattern of magnitude eigenvalues. 
Why should that be so? 

Another difficulty is associated with the prediction of Biferale [llj] Jl2 that there 
would be a change in behavior at the special value of e for which the contributions to 
the helicity sum in equation (6"T5) grow toward the high-n end of the inertial range. 
The value is 

e blf := 1 - A" 2 / 3 . (26) 

We have just seen that the model is stable for e = 0.3, which lies below the Biferale 
value e approximately 0.37. Look at figure |2|, which is the analog of Figure [l| but 
now for e = 0.5. The magnitude spectrum looks much the same as before, but there 
is a qualitative change in the phase spectrum. Now, there are a large number of 
unstable eigenvalues in the phase spectrum. We interpret what we are seeing by 
saying that there is a qualitative change in properties of the asymptotic [y goes 
to zero) model at e = e^/- For < e < euf there are at most a finite number of 
unstable eigenvalues. At e just above euf the system acquires an infinite number of 
unstable eigenvalues. Figure |3| shows the behavior of the spectrum at the Bifarale 
point. Once again the phase spectrum has a qualitatively new character, while the 
magnitude spectrum remains much the same as before. We shall want to understand 
better how this behavior arises from the linear stability analysis. 



3.3 Eigenvalues in the Magnitude sector 

The scaling behavior of the eigenvalues in the magnitude sector is given by the very 
same ideas which we already used in our analysis of u n . As we shall see, some new 
ideas will be required for the phase sector. For this reason, we shall dispose of the 
magnitudes here and move on to a more extended discussion of phase eigenvalues 
in the next chapter. 

The general structure of the eigenstates is illustrated in figure || Here we look 
at right eigenvectors with j = 26. The eigenvalue equation in the magnitude sector 
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IS 



5<$> J 

w n 



Y,(D nm + C nm )5& m . 



(27) 



rn 



For small n the eigenvector is very small. In fact, the eigenvector must decrease 
rapidly as n decreases to enable the eigenvalue term from swamping the right hand 
side of the equation. This decrease is a falloff from a plateau which occurs at 
n about equal to j. At this point, the eigenvalue term and the cascade term in 
equation (|27|) are about equal in size. As n increases still further, the three parts of 
the cascade term each become much larger than the eigenvalue term. To ensure the 
cancellation of the different parts of the cascade term the eigenvector settles down 
to an oscillation with period three, which holds throughout the entire ISR. Finally 
as n enters the VSR, the eigenstate once again falls off quite rapidly. The behavior 
of the phase eigenstate is much the same, except that there is no falloff in the VSR. 

One can describe the same process in physical terms by saying that the eigenvalue 
term adds energy to the system which then cascades toward the VSR, where the 
energy is dissipated. 

Now one can see why the dissipation plays such a large role in determining the 
behavior of the eigenvalues. The deviation 5Q 3 n produces a change in the conserved 
quantities, both the helicity and the energy. This change cascades down through 
the ISR toward the VSR. Because of the conservation of energy this cascade is not 
damped and remains constant in size until the VSR is reached. Since the eigenvector 
remains constant into the VSR, then the VSR behavior serves as a kind of large n 
boundary condition on the eigenstate. 

A simple count shows how this works. For large n in the ISR, <5$ ra is periodic 
with period three. Then, there are two parameters which describe the wave function: 



In our previous work |§ we found that there were two conditions upon the large-n 
velocities, required to keep these velocities from blowing up deep into the viscous 
subrange. Two parameters, two conditions. Everything is determined. Thus we 
might expect that all eigenstates with sufficiently small j would have the very same 
values of the parameters for large enough n in the ISR. Table 1 serves to check this 
point. In this table we have shown values of the ratio of these parameters for various 
different eigenstates. According to the theory, the ratio should be unity. Clearly 
the theory works well for the magnitude eigenstates. (The table also shows that the 
constancy of the parameters does not work for the phase eigenstates, but that story 
will be told later.) 

Thus, the matching into the VSR uses two of the four boundary conditions on 
our linear chain problem. A very similar mechanism sets the other two boundary 
conditions in the regions in which n is of order j. There are two quantities which can 



jp 2 {k) 



Pi 



<5$3fc+2 



(28b) 
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determine the behavior of the eigenfunction in this region. The first is the variation 
in helicity coefficient B produced by the disturbance. This variation produces a 
term in which varies as (e — l) n and which then grows to be of relative order 
unity when n is of order j. By setting this coefficient, and also the value of the a\ 
one has two coefficients at ones command. These two are just enough to ensure that 
the eigenfunction decays, rather than grows, as n goes to one. 

There is only one important difference between the problem of determining the 
static solution u n and the eigenf unctions. For the former, we have forcing on the 
first shell. In the latter, the important forcing is produced by the eigenvalue term, 
and occurs on the j th shell. The scaling analysis for the j th eigenvalue is modified 
because the range of the cascade becomes j to rather than the 1 to No that we 
used in the analysis of the velocity. Thus the scaling rule, with corrections, becomes 
on the small- j end 

a 'U) = a (3)(i + o(i) 4 / 3 )) for j in region S (29a) 

Here, as before, the prime indicates a decreasing in the viscosity by a factor of A 4 
together with a shift in the cutoff, N. The corresponding result on the large- j end 
is the statement 

a 'U+3) = AV j) (l + 0(e - for j in region L (29b) 

Both together give the scaling law: 



a 



(i+3) 



X 2 a^' for j in region I (30) 



To check our analysis, we should check these rules. The first check, done in figure 
0a, is to plot overlays of the primed and unprimed eigenvalues on polar plots like 
that in figure |l[ We show the phase eigenvalues here, but the scaling fits are even 
better for the magnitude eigenvalues. The two spectra overlay precisely at small j, 
as expected, and fit badly for large j. In the second check, one plots overlays of a 
and cr'/A 2 as shown in figure 0b. This overlay shows agreement between the two for 
large eigenvalues but not for small. A more careful check is to take the ratio of the 
nearby eigenvalues in these two figures. Call this ratio R. The quantity \R\ — 1 is 
a quantitative measure of the errors in our statements ( |29a| ) and (|29 b| ) . Fig || plots 



this ratio versus j and shows that the order errors vary as stated in these equations. 
Thus, one can feel that the magnitude response is understood reasonably well. 



3.4 Eigenvalues and Conservation Laws 

One basic principle about this model is that the dissipation can play a large role in 
determining the ISR behavior. To see this fact in more detail, we examine the effect 
of the conservation laws for energy and helicity upon the eigenvalue analysis. 
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Let us recollect that the conservation laws for the the system can be expressed 

as 

J2U n C n {U)h n = 0. (31) 

n 

Here h is a quantity which defines the two conservations. For energy conservation 
h — 1; for helicity 

h={e-iy l . (32) 

For the j th right eigenvector and its eigenvalue <7 J , we have an eigenvalue equa- 
tion of the form 

<T j 5& n = £(±C nm - D nm )5& m . (33) 

m 

Here the minus sign corresponds to the magnitude eigenvalues and the plus to the 
phase eigenvalues. Multiply by (u n ) 2 h n and sum over all n to obtain 

a^M 2 h n 5^ n = ±Y / h n (u n ) 2 C nm (u)5& m - h n (u n ) 2 k 2 n 5& n . (34) 

n nm n 

The conservation identity ([H]) is true for any U and hence also for U = u(l + 5$). 
Since is a small perturbation we can expand around u. 

= ]T h n U n C n (U) = h n [u n C n (u)[l + 6& n ] + J2M 2 h n C nm ( u )6& m }. (35) 

n n nm 

In this way we have expressed the Jacobian ( |19b| ) in terms of the cascade itself. 
(HH, and |) yield 

- h n (u n ) 2 C nm 6& m = h n u n CJ& n = Y,(~Dn + F n )u n h n 5¥ n . (36) 

nm n n 

Now we can rearrange the eigenvalue equation. In the case of magnitude response 
the dissipation term in equation (^) adds to the dissipation term in equation ( j33|) 
to give us an identity for the eigenvalue: 

^E ft "W 2 ^ = -2vY,h n {u n k n ) 2 8& n + fh Ul 5&i. (37a) 

n 

The left hand side of this equation is the rate of decay of the conserved quantity, 
as determined by the eigenvalue. On the right we see that the decay is (naturally 
enough) not produced by the cascade but only by the dissipation through viscous 
damping and also by the addition through the external force, /. Thus we understand 
once more that the dissipation must have a crucial role to play. 

If we go through the same calculation for the phase response, the result is totally 
different. Instead of adding to one another, the dissipation terms cancel out (!), 
leaving us with the identity 



J2 h n (u n ) 2 5& n = -fh Ul 5${. (37b) 
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Compared to the magnitude sector, the phase sector shows a quite different form 
for the conservation law identities. We can no longer say that dissipation produces 
decay. Instead we say that the relevant quantity is added through the force term 
and then changes in time because of this addition. We now turn to a more detailed 
consideration of the phase sector. 



4 Phase Response 

4.1 Establishment of phase transitions 

From what we have seen, both u n and the entire magnitude sector of the linear 
response vary smoothly as e passes through the Biferale value given by equation 
(ffif). When A = 2, this transitional value is 0.37. The story is different for the 
phase response. Whenever e passes through the critical value, then there is a quite 
apparent change in the structure of the eigenvalues. As the viscosity goes to zero, 
this change involves having a large number of eigenvalues pass from being stable to 
being unstable. In the asymptotic limit, one third of the entire set of ISR eigenvalues 
undergo such a passage at this point. 

To see the evidence for this proposition, return to part b of figures [l],|2|, and figure 
|3|. These pictures respectively apply below, above, and at the phase transition in 
e. Away from the phase transition, the phase eigenvalues fall into two classes: The 
eigenvalues fall onto three lines of constant phase: one real and two with opposite 
phases. Others, the "deviating eigenvalues" do not seem to fall into the simple 
pattern, and have phases which change from eigenvalue to eigenvalue. Since we are 
interested in scaling properties, we want to know something about the pattern which 
arises when the viscosity is taken to zero. An examination of the spectrum shows 
that its three main branches grow longer as v becomes smaller, but that the number 
of deviating eigenvalues does not grow. Figure § presents a counting of eigenvalues 
in the phase response sector. They are divided into the following categories: real, 
"constant phase" complex, and deviating eigenvalues. They are then counted at 
different values of Nd. We see that the number of deviating eigenvalues remains 
the same as we change the number of shells in the ISR. The same analysis -and 
result- applies above the transition. Consequently, for v — > the number of deviant 
eigenvalues becomes negligible compared to the number in the three main branches. 
In the asymptotic limit, each of the main branches has one third of the total eigen- 
values. Thus, the deviating eigenvalues are not a scaling limit phenomenon, but 
rather a transient which defines the approach to scaling at the ends of the ISR. 

There is a main branch on the real axis both above and below the transition. 
Above the transition, this branch is unstable; below it is stable. At the critical point, 
we have a change involving, in the v — > limit, an infinite number of eigenvalues. 

Figure [10], together with figures |l|b,|]b, and presents the flow of phase eigenval- 



ues as e goes from below to above the transition. We start at |T0] a with the familiar 
spectrum at e = 0.2. Real eigenvalues collide and turn to complex ones and these 



13 



(so created) deviating group of eigenvalues wanders towards the imaginary axis. If 
we had chosen a larger system the number of deviating eigenvalues would still be 
the same, and only the straight branches would gain in members. The smaller the 
system is the more is the discontinuous transition smeared out by finite size effects. 
(For finite v also the critical point e c weakly depends on v and slightly deviates 
from €bif, see figure 0. These finite size effects spoil the regularity of the phase 
transition. A related transition with less disturbing finite size effects is discussed in 
appendix B) 

As the deviating group is close to the phase transition the two conjugate straight 
branches of eigenvalues split into two (Fig. causing a breakdown of scaling law 
S. A scaling of the form S 2 (i.e., make use of the symmetry n — > n + 6) however, 
is still valid. The scaling of the solution and of the magnitude eigenvalues does not 
break down at this point. 

Above the transition fig |10| d the deviating group has curved in the other direction 
and the eigenvalues return to the now unstable side of the real axis. 

There is a hidden structure of the deviating group not easily seen in our repre- 
sentation (^); these eigenvalues fall on a straight line in a log(\a\) — arg(a) plot. 
Figure |11] shows (the upper half) of the deviating group and the real eigenvalues for 
e = 0.33, 0.37, and 0.4. A dotted line marks the border of instability at it/2. Remem- 
ber again that only branches with constant phase (here horizontal) will contribute 
in the iV — > oo limit. From figure [11] it is evident that with N — > oo instability 
will occur immediately above e^f, since an arbitrary small tilting suffices to produce 
unstable real eigenvalues. Moreover, for v — > there is not even a finite number of 
unstable eigenvalues for e < euf (see Fig. [12]). This means that euf coalesces with 



tbif — e c - (3£ 



Although the instability mechanisms consists of one pair of complex conjugate 
eigenvalues after another crossing the imaginary axis, at iV — > oo an infinite number 
of such pairs coalesce and an infinite number of oscillatory instabilities is unleashed 
within an arbitrarily small change in e. 

For v — f the destabilization scenario of the fixed point is thus different than 
has been suggested in the light of "finite-size simulations" (|J [|). On one side of 
the phase transition we have totally stable behavior, while on the other side we have 
immediately an infinite number of unstable modes on all scales. These instabilities 
are purely exponential and not oscillatory. 



4.2 Asymptotics 

The phase transition is a change in the behavior of a very large number of eigen- 
functions all at once. Biferale JTTJ has explained this phase transition as a blockage 



in the energy flow caused by a flow of helicity. Since the blockage can occur any- 
where in the ISR it seems reasonable to assume that, when the conditions are right, 
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blockages can occur at many places affecting many eigenf unctions and eigenvalues 
at once. 

Another way of asking the same question is more mathematical in character. To 
see this phase transition we must have some kind of change in the small-i/ asymp- 
totics of the system. That is we must have some sort of change which can be seen 
by an eigenfunction with j much larger than unity and much smaller than Np. So 
we need an asymptotic theory of eigenfunctions for this system. 

To understand the phase transition, we must first understand how the phase 
sector can be so unstable. The key is given in Table 1. In that table we see that 
in contrast to the magnitude sector, the phase sector shows far more flexibility in 
the large-n limit of the eigenfunction. The magnitude sector has but one value of 
Pi and P2 for all eigenfunctions in this region. The phase sector can actually have 
two different linear combinations. The two permitted combinations are: 



5$ 3fc+ i = <J$ 3 fc+2 = -5$3* (39) 



-5® 3k+1 = 5<5> 3k+2 = <5$ 3fc (40) 



and also, for example, 



for integer values of k. The changes in expression (|39|) are an exact symmetry of the 
GOY system and produce a phase sector eigenfunction which has eigenvalue zero. 
The changes in expression ( fiO"D satisfy the eigenvalue equation for all n-values except 
n = 1. Here, the eigenvalue equation fails because of the forcing term. However, 
this combination also forms a possible high-n behavior. All eigen functions have one 
or another linear combination of these two behaviors as the high-n limit. But there 
are two such linear combinations possible. In comparison to the magnitude case, we 
seem to have lost a boundary condition at the high-n end. 

Thus, we state the first contrast between the magnitude sector and the phase 
sector. The magnitude sector can be described by giving two boundary conditions 
at the high-n side of the ISR; the phase sector can be described by giving only one. 

We must have someplace an extra boundary condition for the phase sector. Bi- 
farale directed our attention to the conservation of helicity. Instead of doing a full 
analysis of the possible extra conditions, we just follow his direction and look at 
the helicity conservation law, equation (|37b|) , in the phase sector. (The reader will 
recognize that some leap of faith is required around this point in the argument.) 
This equation is 

I»„m; = -^2i. (4i) 

But, for an eigenvalue with j in the middle of the ISR, the first component of the 
eigenfunction is very, very small. This component goes to zero more rapidly than 
an exponential in j. For this reason, it is quite reasonable to neglect the right hand 
side of the helicity identity and write instead: 

h n (u n ) 2 5& n = for j in region I (42) 
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We propose that this identity replaces the lost high-n boundary condition for the 
asymptotics of the phase sector. 

This proposition has several quite attractive features. We expect the eigenfunc- 
tion to be of roughly the same order of magnitude over the whole range in which n 
is greater than j. Then for e < e^f the summation has its main contribution at the 
lower-n end and then falls off as 

8& n (-l^f) n . (43a) 

When we are at the critical point, the summation does not fall at all but the sum- 
mand looks like 

M£(-1) B . (43b) 

Finally, above the critical point the major contributions come at the high-n end and 
then fall off toward lower n with the same law as shown in equation (|43a]) . Near 



the phase transition, these behaviors produce long-range correlations between the 
different parts of the ISR. It is these correlations which form the key to the phase 
transition. 

From these assumptions, one can find scaling laws for the phase eigenvalues, 
namely on the small- j end 



a' 3 = a j 



/l_ e .\JVb-i\ 
1 + f — — J j for j in region S (44a) 



Here, as before, the prime indicates a decrease in the viscosity by a factor of A 4 
together with a shift in the cutoff, N. The corresponding result on the large- j end 
is the statement, which we have used before. 

0./G+3) = X 2 a {j \l + 0(e - for j in region L (44b) 

Both together give the scaling law: 

a = X 2 a^ for j in region I (45) 

To check the accuracy of our thinking we show in figure |7| polar plots in which 
we overlay the eigenvalues <jj with, respectively, the eigenvalues cr'- and A~ 2 aj +3 . 
According to the theory, the first of these should agree very well for small magnitudes 
of the eigenvalue, while the second should agree on the large magnitude end. The 
figures bear this out. For a more accurate check, we construct errors as, for example, 



5 



1± 



1. (46) 



These errors are shown in figure || and compared with the theoretical estimates taken 
from equations ( 44b|) and ( f44a| ) . Since the estimated and the actual errors agree 



rather well, we can argue that we have caught the essence of the phase eigenvalues. 
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5 Summary and Conclusions 



The GOY model has a linear response behavior about a static solution in which the 
stability eigenvalues extend over a huge range of frequencies. In the ISR, the eigen- 
values show a scaling behavior limited by disturbances from stirring or from viscous 
effects. The scaling of the eigenvalues is more complex in the phase sector than 
in the magnitude response. The phase sector shows a phase transition connected 
with a boundary condition which moves from one end of the ISR to the other. This 
boundary condition is derived from the conservation law for the model's version of 
helicity. Eigenvalues change quite abruptly as the model's parameter passes through 
its phase-transition value. 

Much of the scaling mirrors the scaling properties of the static GOY model. 
However, there is scaling associated with the phase transition which we have not 
investigated in any detail. We have also not fully established the scaling behavior 
of the linear-stability eigenvectors. These studies are left for the future. 

For now, we have demonstrated scaling in the linear response of a still system. 
We have seen a new phase transition and gained a qualitative understanding of its 
source. 
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A Large scale forcing 



Most studies of the GOY model || [| [| employed F n = f5 n ^ as large scale forcing. 
Though this forcing seems innocent, it has some disadvantages. The velocity com- 
ponent U3 is much smaller than 1*1,2,4,5,6,... Hi? what is clearly unphysical. We force 
the system on level n = 1. As can be seen from equations (0) and (|), with w 3 = 
and n — > n — 3, the shift in the forcing by three shells just shifts the solution by 
three shells as well. Q Q 

With the different forcing the borderline of stability becomes only slightly dif- 
ferent. For instance, e c beyond which (^) becomes unstable is in the standard case 
(A = 2, v = KT 7 , / = V2 ■ 5 • 10~ 3 ) shifted from 0.370 to 0.377. More significant 
is the change in the amplitude of the e c oscillations with viscosity. This viscosity 
dependence is shown in Fig. [12] for A = 2, and looks similar for other values of A. 
In contrast to the case with forcing on the 4th shell (see |§, Fig. 12) the stability 
border e c now seems to approach the constant e^/ for vanishing viscosity v — > 0, 

lime c = e fci/ := 1 - A~ 2/3 . (47) 



Some of the transitions to instability in Fig. [12] are due to a complex pair of eigen- 
values passing through the imaginary axis, others are due to a single real eigenvalue 
going through zero. f\ 

In the main text we have studied the behavior of the overall eigenvalue spectrum. 
The stability border, on the other hand, is determined only by the first eigenvalue 
crossing the imaginary axis. Nevertheless, in the main text we have two strong 
evidences for a dependence of the stability border on the particular forcing. As we 
have seen in section 3.2 the modulus of the eigenvalues is set by the cascade term. 
The eigenvalues that decide the stability are the small ones close to zero. Their 
size should depend on how the solution looks like at the very first shells, and this 
is obviously influenced by the kind of forcing. Also formula ( |37b| ) shows that the 
stability border should depend on the particular forcing. 



4 With the forcing F n oc 8 n .i a static solution of eq. (g) exists for all < e < 1, i.e., the 
saddle node bifurcation found in || is an artifact of the forcing F n oc <5 n ,4- In particular, the 
overall existence of the static solution now allows for an analysis of the eigenvalue spectrum for 
the parameter values X = 2 and e = 0.5 [|] |] f|. 

5 The artificially small u 3 (for F n oc <5 n! 4-forcing) also leads to strongly non-normal eigenvectors. 
This was first discovered by J. P. Brunet JlQj] , who also calculated the corresponding pseudoreso- 



nance. The reason for that can be seen from (19b, 20b): (A,p)3 im with m = 1,2,4,5 is small as 
1*3 ~ v is so small. All the eigenvectors of then have a huge 3 rd component and are thus nearly 
parallel to each other. This makes very nonnormal, the more, the smaller v is. 

6 The model with the forcing on the fourth shell shows a transition from stable to unstable 
behavior only via two complex conjugated eigenvalues going through the imaginary axis. With 
the forcing on the first shell instability for some v occurs via a real eigenvalue going through zero, 
followed by a series of Hopf bifurcations (take e.g. v = 10~ 8 , / = v^-0.005, k = 0.0625). However, 
from the viewpoint of the scenario of instability in the limit v — * 0, treated in section 4.1, this 
additional type of transition is not fundamentally different from the other one. 
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B A related phase transition 



In this appendix we analyze the related phase transition in the spectrum to the 
matrix 

A 5 = £ • C - D. (48) 

The matrices C and D contain the solution of (fj). For general £ the matrix is 
not the stability matrix of the system. Only for £ = — 1 we have A_i = Am and 
for £ = 1 we have A\ = A$. The generalized amplitude matrix A^ with £ < 
shows the same features as A M itself. So we only study the interesting case of A^ 
with £ > which displays a phase transition at £ = 1. The advantage of studying 
this phase transition rather than the one in e discussed in the body of the paper is 
that here the critical point £ = 1 does not depend on the viscosity. Thus finite size 
effects are less complicated than for the transition in the spectrum of A$ at the v 
dependent e c . Taking £ ~ 1 (but £ ^ 1) can be thought of as equivalent to slightly 
varying the dynamical equation (0) and thus changing the phase symmetry of that 
equation, which leads to a different solution. We perform the analysis for e = 0.3 as 
it is the spurious stabilization of the phase matrix for e < e c « 0.37 which we want 
to understand. 

We start with £ = 0. All eigenvalues — vk\ of A = —D are of course stable. 
However, with increasing £, more and more eigenvalues of A^ turn positive as soon as 
£ > v. This feature holds for | lgz/| ~ 30 orders of magnitude in £. From comparing 
the viscous and the inertial contribution to A^ we obtain that the largest eigenvalue 
increases as 

a max ~ e ,2 v- 1/2 - (49) 

o~max becomes as large as 10 12 for £ « 0.63 and the standard parameters A = 2, 
e = 0.3, iV = 81, v ~ 10~ 30 . This is the behavior we see in the spectrum of A^ = Ax 
for e > e c and what we had expected also for smaller e. 

Why and how is stability achieved for £ = 1, i.e., why does A\ = A^ not have 
any positive eigenvalue for e < e c ? 

For £ growing beyond 0.63 towards 1, two real eigenvalues merge, form a complex 
pair which moves in the complex plane on a circle and finally turns stable through 
an inverse Hopf bifurcation. This happens again and again through a selfsimilar 
cascade of bifurcations towards phase symmetry at £ = 1 - no unstable eigenvalue 
is left. This phase symmetry of A\ = A^, discovered in |J, reflects itself in a zero 
eigenvalue, as discussed above. 

We now analyze the singularity in detail. We introduce the distance r from the 
singularity at £ = 1, 

£=l + r (50) 

and increase |r| on a log scale from r = 0. In this way we break the phase symmetry 
in a stronger and stronger way. 

We describe the features for r < 0. For r = we only have negative eigenvalues 
and one eigenvalue equals zero. This center manifold eigenvalue turns positive for 
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\r\ > 0, signaling instability. Moreover, two different small (modulus wise) real 
eigenvalues (< 0) merge on the real axis and form a complex pair which wanders 
towards the 9?er = axis and turns unstable via a Hopf bifurcation. This happens as 
early as for r m 10~ 10 , see figure O. In the right complex half plane they continue 
to wander on a kind of circle, finally merging again on the real axis, now forming 
two positive real eigenvalues. But meanwhile a second pair of complex eigenvalues 
has formed in the left half plane which again turns unstable via a Hopf bifurcation 
and then finally merges. This happens again and again, leading to 3, 5, 7, 9, 11, . . . 
positive eigenvalues. The real parts of the positive eigenvalues are plotted in fig. 



13| . The selfsimilarity of the cascade of bifurcations can clearly be observed. The 
maximal eigenvalue a max oc r 4//3 . The scaling law mirrors the classical K41 scaling 
of u n . 

For r > the behavior is very similar to the r < case. The only difference 
is that the center manifold (zero) eigenvalue turns stable first and we now have an 
even number 2, 4, 6, 8, . . . of positive eigenvalues. 

To demonstrate the finite size effects of the transition we plot the number of 
positive eigenvalues as a function of lgr for various system sizes N, i.e., viscosities 
v, see figure In the iV — > oo limit, even for an arbitrarily small \t\, we have an 
infinite number of unstable eigenvalues. 
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Table 



n 



magnitude 



phase 



Re( 



(n)< 



7+T 
Pi 



Re{ 



p 2 ( n ) ' 

pj, +1 (n) ' 



p 2 ( n ) \ 



18 
18 
18 
18 
18 
18 



35 
45 
50 
55 
62 
70 



1.02897 
1.00132 
0.99986 
1.00004 
1.00000 
1.00000 



0.99597 
1.00041 
1.00002 
1.00002 
1.00000 
1.00000 



0.02867 
-4.83889 
0.02867 
0.31876 
0.02867 
0.31876 



0.81927 
22.2840 
0.81891 
-0.00242 
0.81891 
-0.00242 



13 
13 
13 
13 
13 



16 
35 
45 
55 
71 



75.5759 
0.99648 
0.99984 
1.00000 
1.00000 



2.08638 
1.00048 
0.99995 
1.00000 
1.00000 



-0.58804 
40.53829 
-0.00677 
-0.42772 
40.53836 



-51.81228 
-0.11043 
0.02073 
-51.26321 
-0.11044 



Table 1: Behavior of parameters for eigenstates. Here j refers to the j th eigen- 
function while n describes the components of that eigenfunction. These two pa- 
rameters describe the large-n behavior in the ISR. They are fixed by the value of 
v = 10 -3 16 -21 in the magnitude sector but vary considerably in the phase sector. 
The ratios are essentially the same for the n-values in between the given values. 
(Correlations between phase values with n different by a multiple of 3 stem from 
the period 3 in the eigenvectors.) 
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Figures 



Figure 1: Linear Stability Eigenvalues. The case considered is e = 0.3, v = 
10 _3 16 -26 , N = 90. Here we plot the amplitude (diamonds) and phase (circles) 
of the eigenvalues in a kind of polar plot. The phase in this polar diagram IS the 
phase of the eigenvalue. The radial coordinate is given by equation ( p5[) which 
essentially produces the logarithm of the eigenvalue. 



Figure 2: The same as Figure [I] except that now e has the value 0.5. The amplitude 
matrix becomes eventually unstable at e « 0.558. 



Figure 3: The same as Figure |l]b except that now e has its critical value, euf which 
is about 0.37. Notice that there are now six branches, all with complex eigenvalues. 
The equal spacing on each branch is evidence of scaling. But, evidently, the scaling 
is quite different from that shown in figures |I|b and |2|b 



Figure 4: We plot W n = k]l z V n against n. Notice how W n oscillates in the inertial 
range and how it falls off quite rapidly in the dissipative range. (The latter is for 
n > Nd « 67.) The actual calculation is done for e = 0.3, v = 10~ 3 16~ 20 (circles) 
and v = 10~ 3 16~ 21 (dots). The behavior of W in these two simulations is compared. 
In the second (dots), N is increased by three and v is decreased by a factor of A -4 . 
For all n smaller than 56 or so the two calculations agree, (b) A comparison like 
that in a) except that now W n is compared with W^_ 3 . These agree top plotting 
accuracy for all n bigger than 20 or thereabouts. 



Figure 5: The error in the scaling relations for the velocity. The deviation from 
unity of the ratio of two pieces of data in equations (|16a|) (triangles) and ( |16b|) 
(diamonds). The slope of the solid lines show the theoretical estimates of the error. 
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Figure 6: The general structure of right eigenvectors for the the amplitude 
5^u n k^ 3 (solid) and phase 5<S>^ (dashed line) (N = 75, e = 0.3, v = 10- 3 16- 21 ). 
In both cases the modulus of the eigenvector is displayed. The amplitude eigenvec- 
tor gets damped when it enters the viscous range, while the period 3 of the phase 
eigenvector continues into the viscous range. 



Figure 7: a) Eigenvalues of the phase matrix for e = 0.3, N = 90, and N = 93. One 
can see three main branches and the validity of scaling law S. b) The same data, 
but the eigenvalues for N = 93 are shifted by 2. Scaling law L is at work. 



Figure 8: Errors in the scaling laws for the eigenvalues, equations ( |29aj , |44a| ) and 
( [29b , 44b ). Here we plot \R\ — 1 where R is the ratio of the two sides of the equation. 
The diamonds describe the ratio in equations ( |29a| , [44a| ) while the triangles do the 
same for equation (|29b| , |44b | ) . The unfilled symbols are for the amplitude matrix 
and the filled ones stand for the phase matrix. Lines are theoretical estimates of 
the errors. Part (a) is for e = 0.3 and part (b) for e = 0.5. The errors are for the 
modulus of the eigenvalues, but the errors in arg(R) follow the same trends. 



Figure 9: Number of real (solid line), "straight" complex (dashed), and deviating 
eigenvalues (dots) with increasing shell number The slope of the lines is exactly 
2/3, 1/3 and respectively. 
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Figure 10: This sequence of pictures shows the phase eigenvalues with changing 
e. We start out with three completely stable branches in (a). Real eigenvalues 
turn complex and move towards the imaginary axis (b) until fewer and fewer real 
eigenvalues are left (c). Notice, all the 'deviating' eigenvalues that do not lie on 
a straight line are finite size effects, i.e. they constitute only a negligible part 
for a sufficiently large matrix. These deviating eigenvalues form a straight line at 
e bif — 0.37 (Fig. |3|). Now the spectrum consists of six complex branches and the 
scaling law (S) is no longer valid and has to be replaced by scaling law (S 2 ). Above 
ebif (d) eigenvalues return to the real axis, but are now unstable. Essentially there 
are two complex branches and one real branch again as it was below e^/. An even 
higher value of e = 0.5 is shown in figure 0b. 



Figure 11: Main branch of 'deviating' eigenvalues and real eigenvalues in a log\\a \ — 
arg(cr) plot, for e = 0.33 (triangles), e = 0.37 (squares), and e = 0.4 (circles). The 
tilted branch rotates with changing e and is horizontal at the phase transition and 
unstable above (black). It can be seen that in the N — > oo limit an infinite number 
of eigenvalues turns unstable immediately above the phase transition. 



Figure 12: Border of stability e c as a function of v for A = 2. e c approaches the 
Biferale prediction (thick horizontal line) for small v. Some of the points are due to 
a complex pair of eigenvalues passing through the imaginary axis, others are due to 
a single real eigenvalue going through zero. 



Figure 13: All positive eigenvalues of for r = £ — 1 between and —1 on a 
log- log scale in |r|, corresponding to a lg vs lglg£ plot. The selfsimilar cascade 
of bifurcations towards phase symmetry has its origin in the self similarity of the 
matrix itself. The parameters are N = 80, v ~ 2 • 10~ 29 , e = 0.3, A = 2. For r = — 1 
(right edge of the plot) we have = —D and all eigenvalues are again stable. 
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Figure 14: Number of positive eigenvalues of Ai_ T as a function of \r\ for various 
system sizes N = 80,68,56,44,32,20, left to right, corresponding to viscosities 
v — u ■ 16*, % — 0, 4, 8, 12, 16, 20, respectively. u ~ 2 ■ 10~ 29 , e = 0.3, as in figure [13. 



*e-mail: lohse@cs.uchicago.edu 
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